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Abstract 

An efficient matrix formalism for finding power series solutions 
to boundary value problems typical for technological plasticity is de- 
veloped. Hyperbolic system of two first order quasilinear PDEs that 
models two-dimensional plastic flow of von Mises material is converted 
to the telegraph equation by the hodograph transformation. Solutions 
to the boundary value problems are found in terms of hypergeomet- 
ric functions. Convergence issue is also addressed. The method is 
illustrated by two test problems of metal forming. 

1 Introduction 

A model of incompressible, ideal, rigid-plastic material was developed about 
a hundred years ago. In the two-dimensional case the model is described by 
a hyperbolic system of quasilinear equations. In 1920s Prandl, Hencky, and 
Mises, among others, suggested a systematic way, named slip-line method 
(SLM), for finding stress fields and associated velocity fields for two-dimensional 
plain strain deformation problems for rigid-plastic body . The method is pre- 
sented in many classic books on plasticity theory (cf. pQ, [2], [4]) along with 
its applications to metal forming processes and problems of plastic failure 
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(limit analysis). Solution of a 2D rigid-plastic flow problem by this method 
is reduced to construction of the field of characteristics (slip-lines). Assumed 
domain of the plastic flow is to be decomposed into a set of subdomains 
(patches) of a priori specified types, and appropriate boundary value prob- 
lems (BVPs), consistent with physics of the process, are stated for each of 
the subdomains. The subproblems were usually solved numerically. In 1980s 
a matrix implementation of the SLM, called matrix method, was developed 
within engineering community [5]. The matrix method is based on the as- 
sumption that boundary data are defined by real analytic functions. The 
solution to each subproblem is obtained by applying appropriate matrix op- 
erators to the vector(s) of coefficients of the data. In [6] a new, simpler 
approach to implementation of the matrix method was developed, mathe- 
matically justified, and equipped with recurrent definitions of five matrix 
operators needed for solving main boundary value problems typical in tech- 
nological applications (versus twenty operators in the original engineering 
version of the method). It was also shown that solutions to the subproblems 
can be written in terms of hypergeometric functions. 

In this work explicit expressions for the five matrix operators of the slip- 
line analytic technique (SLAT) developed in [6] are presented. Convergence 
issues are addressed and model examples are provided. 

SLAT can be used as a source of test problems for numerical methods 
aimed at problems with discontinuities. The fact that solutions to problems 
in plasticity may have singularities along characteristics is well known and 
experimentally confirmed. Presence of discontinuities is an important feature 
that cannot be easily detected by numerical methods. In general, capturing 
singularities is very challenging mathematical problem that is of great in- 
terest in many applications, including phase transitions and microstructure 
formation. SLAT also provides examples of analytically solved problems with 
free boundaries since only a structure of the plastic domain is "guessed" a 
priory and finding the boundary of the domain is part of the solution. Fur- 
thermore, although analytic solutions have been superseded with powerful 
numerical methods, in particular FEM, the slip-line analysis is a working 
tool in metal forming [5] , granular flow modelling [3] , and geomechanics [TD] . 
The velocity field in the hodograph plane can also be found by this method 
[6]. Exact solutions for stress and velocity fields can be used to derive ana- 
lytic expressions for the plastic power, which provides a framework for finding 
optimal in a certain sense geometric or material parameters of the process 
under consideration. 
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SLM is also a classic topic in standard courses on plasticity theory, and 
so is the method of characteristics in PDE theory. Therefore the authors 
believe that an elegant matrix formulation of the method presented in this 
paper has a good potential in engineering and mathematical education due 
to simplicity of the material model, lucidity of the basic ideas, and classical 
character of mathematical machinery. 

The paper is organized as follows. 

In Section 2 the governing equations for von Mises ideal rigid-plastic mate- 
rial under conditions of plain strain deformation are introduced and a classic 
reduction to quasilinear hyperbolic system of PDEs is given. Using hodograph 
transformation, the system is transformed into the telegraph equation, which 
means that the original system is C-integrable. This fact was known long ago 
[2], but here the transformation to the linear problem is done by appropriate 
change of variables and does not involve any considerations from mechan- 
ics. In Section 3 three main BVPs typical for applied plasticity are defined. 
Using Riemann function for the telegraph equation, the initial characteristic 
problem with analytic data can be solved in terms of Bessel functions [8]. 
We show that this solution can be expressed in terms of the hypergeometric 
functions depending on the product of independent variables, which is more 
practical for computations. Exact solutions for Cauchy and mixed problems 
are found by reducing the problems to equivalent initial characteristic prob- 
lem. In Section 4 the method is applied to two test problems from the theory 
of plastic flow: 

• Calculation of the slip-line field generated by two circular arcs, and 

• Computation of the stress state near an elliptic hole loaded with con- 
stant normal pressure. 

Solutions obtained for the test problems are in excellent agreement with those 
known from literature. 

To simplify exposition, the following notations and conventions will be 
used hereafter. 

• t n = t n /n\, where t is a variable. 

• [t] = [l,ti,t 2 ,t 3 , . . .} T (vector-column). 

• Any given function is real analytic (i.e., it can be represented as the 
sum of a power series). 
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Summation is taken over the range from to oo unless specified other- 
wise. 



• Notation 

oF 1 (n + l;z)=n\J2 k]{k + n)] 

stands for the hypergeometric functions. The subindexes will be omit- 
ted for simplicity. 

• Given a power series ^ a,k tk, notation a stands for the row of its coef- 
ficients. 

2 Governing equations 

Let a homogeneous isotropic body made of ideal rigid-plastic material be in 
a state of plane strain plastic deformation. Assume that the flow is parallel 
to rry-plane. Let D denote the projection of the body onto this plane. It 
is assumed that the stress tensor components o"^, i,j = 1,2, and the strain 
rate tensor components e^- in the domain D satisfy the following equations: 

1. Equilibrium equations 

a^ = 0, J = 1,2. (1) 
Here the Einstein rule of summation over repeated indices is assumed. 

2. von Mises yield criterion 

{^-^f + Aa\ 2 = Ak\ (2) 
where k is the shear yield stress of the material. 

3. Constitutive relations 

= , J — 1,2. (3) 

012 e 12 

4. Incompressibility constraint 

en + e 22 = 0. (4) 
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Two perpendicular lines, called slip-lines, are passing through any interior 
point of the plastic domain. Each of these lines is tangent to a maximum 
shear stress direction at this point. Lines in the first and the second shear 
directions [7] are called a- and /3-lines, respectively. Under the plane strain 
conditions, the mean stress a is given by a = (an + 0-22) The stress 
components at a point (x, y) are determined by the equations 

cr n = a — k sin 20, cr 2 2 = °~ + k sin 20, 0^2 = k cos 2(f), (5) 

where = 0(x, y) is the angle from the x-axis to the a-line passing through 
the point. 

Substituting these equations into ([I]) yields 

a x — 2k (<f) x cos 2cj) + (fi y sin 20) = 0, (6) 
o y — 2k (4> x sin 20 — <p y cos 20) = 0, (7) 

where lower indices indicate corresponding partial derivative. Quasilinear 
system (3) is hyperbolic, and its characteristics coincide with slip-lines. It 
is shown below that the system can be converted into a linear PDE by a 
hodograph transform reversing the roles of the dependent, (cr, 0), and inde- 
pendent, (x, y), variables. 

Assuming that the Jacobian 

J = a x (j)y — (Ty(j)y 

does not vanish, one derives 

4> y = J x a , X = - J y (j , o- x = Jy<t> o-y = - J x^, 

and the system ([7]) transforms into 

x$ = 2k (x /j cos 20 — )— y (j sin 20) , 
y^ = 2k (x a sin 20 - y a cos 20) . 

Introduce point-dependent rectangular coordinate system with axes di- 
rected along characteristics, 

X = x cos + y sin 0, Y = — x sin + y cos 0, (8) 

and the characteristic coordinates (a,/3), 

a = 0/2 + (a - a )/(4k), (3 = 0/2 - (a - a )/(4k), (9) 
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Figure 1: Curvilinear characteristics rectangle with initial a-line OA and 
initial /3-line OB. 



where do is the value of the mean stress at the origin. Then, after algebraic 
simplifications, the governing system of equations takes the form 

Y a + X = 0, X p - Y = 0. (10) 

Geometric meaning of the characteristic coordinates is clear from Fig.l. If 
the functions X, Y are smooth enough, each of them satisfies the telegraphy 
equation 

<Jfe + ' = °- <"> 

The same equation holds for the radii R, S of a- and /?-line curvatures, 
respectively, and for the components of the velocity field [TT] . The stress 
components at a point (x(a, /3), y(a, /3)) are defined by (JSJ), with 

= a + /3, a = a Q + 2k(a- (3). (12) 

The main stresses are given by the equations <j\ = a + k, a 2 = ex — k. 
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3 Main boundary value problems 
3.1 The initial characteristic problem 

It is well known that there exists a unique solution to the following initial 
characteristics problem [7]: 

Find a function f(a,/3) satisfying 07]) in the domain OACB (Fig.l) and 
the initial conditions 

f{a, 0) = ^2 c n a n , a e (0, «i), /(0, p) = d n (3 n , (3 e (0, ft). 

The coefficients c n , d n are given real numbers and the compatibility condition 
Co = do holds. A classic solution in terms of Bessel functions [8] can be 
rewritten in terms of hypergeometric functions as 

f(a, P) = J2 ( c ««« + d M F(n + 1; -a/3) - c F(l; -a(3). (13) 

For a plasticity problem with curvatures of the initial characteristics given 
by the equations 

R{a, 0) = J2 a n®n, S(0, P) = J2 b nPn, (14) 

missing data R(0,j3), S(a,0) are obtained from equations (fTUj) (with X, Y 
replaced with R, S) 

R(0,P) = ^2b n p n+1 + a , S(a,0) = ^a n a n+1 + b . (15) 

Then it follows from ([TBI that the slip-line field in the characteristic rectangle 
OACB is defined by the curvatures 

{a n a n + bn-tPn) F(n + 1; -af3), (16) 
S(a,p) = ^ (-On-iOn + b nPn) F(n + 1; -a/3), (17) 

where a_i = = 6_i = 0. 

The following theorem summarizes these results. 

Theorem 1 Let two arcs of intersecting slip-lines be given by the equations 
p^ ). Then the slip-line field in the characteristic rectangle generated by the 
two slip-lines is uniquely defined by series l[Tb)) , ( f77| ). 
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Furthermore, if the coefficients a^, bk are bounded, the series converge as 
exponential series. If for some q > 



\a n \ < q n n\, 



\b n \<q n n\, 



then for any p, < p < 1/q, the series converge as geometric series with 
ratio r = pq provided \a\ < p, \/3\ < p. 

The estimates of the rate of convergence follow from the inequality 

\F(n + l;z)\ <exp(|z|). 

Remark. A particular slip-line field generated by two circular arcs is often 
used in technological plasticity. In this case, R(a,0) = a , S(0,/3) = b and 
solution takes the form 

R = a F (1; -af3) + b j3 F (2; -a/3), S = b F (1; -a(3) - a a F (2; -af3). 

Solution ffl6l) . ffTTj) to the initial characteristic problem can be written as 
double power series 



R(a,P) 
S(a,f3) 



n.k 



a n o: n+ k(—P)k + b n f3 n+ k+i{— oi)k, 



/ J -a n a n+k +i{-(3)k + b n f3 n+k {-a)k, 



n, k 



or in a matrix form as 

R(a, (3) = (aA(fl) + bB((3)) [a], S(a, (3) 



-a.B(a) + bA(a)) 



(18) 
(19) 

(20) 



Here a, b are the vectors of coefficients in the initial conditions ffl4|) and the 
matrix-functions A, B are defined by the formulas 



A(t) 



B(t) 



1 


-t 


t 2 


-h 





1 


-t 


t 2 








1 


-t 











l 


h 


-h 




-h 


h 


-t 3 
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-h 
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-u 
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-h 


u 


-t 5 


k 


-t 7 
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Figure 2: Characteristic triangle with Cauchy data on non-characteristic 
curve OC. 

If an initial slip-line shrinks to a point, the characteristic rectangle de- 
generates into a fan with slip-line field defined by (1201) with appropriate row 
of coefficients set to zero. 

3.2 The Cauchy problem 

Consider a smooth non-characteristic curve OC. Let the radius of curvature 
of OC be J* (7), 7 being the angle between the curve and x— axis. Given 
two functions (7(7), 0(7), the problem is to find functions a(x, y), <fi(x,y) 
satisfying equations (jSJ), dZJ in the characteristic triangle OAC (Fig. 2) and 
boundary conditions cr = a(j), = 0(7) on OC. There exists a unique 
solution of this problem in the curvilinear triangle OAC [7]. 

This BVP can be reduced to the equivalent initial characteristic problem 
on the fictitious characteristic rectangle OACB. Coefficients a,k, bk for the 
initial slip-lines are found as follows. From an infinitesimal characteristic tri- 
angle with hypotenuse lying on OC (Fig. 2) one has the following elementary 
identities 

ds a J cos(r?) = dspj sin(^) = r(-f)dj, (21) 

where ds a , dsp are the length differentials of the bounding a- and /3-lines, and 
i] = 7 — (p. Using (4), functions a = 0(7) and (3 = f3{pi) can be determined 
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at all points of OC. Substituting the expressions 

ds a = R{a(j),i3(j))da, d Sf3 = -S(a(j), (3(-?))dp (22) 

with R, S defined by ffT6l) . ffTTl) into fl2Tl) . one obtains a system for the un- 
known coefficients (3k- A particular case when both the normal stress a n 
and the tangential stress r n are constant along OC is detailed below. In this 

case 

a = j3, i] = 7r/2 — 0.5 cos _1 (r„/fc) = const, 7 = r\ + 2a. (23) 
Writing the radius of curvature of OC as a function of a, 

r(2a + 77) = r[a], 

one obtains from (T2"T1) 

i?(a, a)/ cosr^ = —S(a, a)/ sin 77 = ±2r[a]. (24) 

These identities (with R, S determined by (|T6l) ). (|T7|) ) imply after some 
mathematical manipulations 

a = 2rC, b = 2rD, (25) 

where the upper triangular matrices C, D, derived from the recurrent rela- 
tions in |6], read 
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For simplicity, notation c, s have been used for cos r) and sin 77, respectively. 
Formulas (I25p determine initial data for the initial characteristic problem 
equivalent to the original Cauchy problem in the domain OAC. 

The solution for Cauchy problem is tested below on a classic example 
with known slip-line field formed by logarithmic spirals. 

Example. Consider a circular arc of radius one loaded with constant 
normal pressure and zero tangent stress component. This implies 77 = 7r/4 
and r = (1, 0, 0, ...). It follows from ([25]) that 

a= V2(l, 1, 1, ...) b = <s/2(l, -1, 1, -1,...). 

Therefore the radii of curvature for the initial a- and f3-Ym.es are R(a, 0) = 
exp(a) and S(0,P) = y/2exp(-0)). These are the radii of curvature of 
logarithmic spirals. 

Stress-free surface boundary (cr n = r n = on OC) can be analyzed 
similarly. However, in this case the shape of the free boundary is not known a 
priori, and to start a process of solution based on the identities ( JUJ) one needs 
to know one of the initial slip-lines. Assuming for certainty that the initial 
slip-line OA (Fig. 2) is known (for example, vector a could be found from 
solution of BVP on the adjacent patch), one finds the row b of coefficients for 
the initial /3-line and free surface radius of curvature r(a) from the identities 
( |24|) . It can be shown that 

b = aF, r{2a + r]) = R(a,a)/V2, (26) 

where the entries of the matrix F are defined as 

/« = -1, f l3 = 0, forj < i, f l3 /2 = far j > i, (27) 

and R(a, (3) is given by ffTB"]) . Finding the shape of a free surface may be of 
interest in applications. 

3.3 Mixed problem 

Consider a smooth non-characteristic curve OC. Let the radius of curvature 
of OC be r(7), 7 being the angle between the curve and a;— axis. Given the 
radius of curvature for the initial a-line OA and a function 0(7) on OC, the 
problem is to find functions a(x, y), 4>( x ^y) satisfying equations (jSJ), (JTj) in 
the characteristic triangle OAC and boundary condition cf) = 0(7) on OC 
(Fig. 3). It is known that the problem has a unique solution [7]. 
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Figure 3: Characteristic triangle for the mixed problem. One of the bounding 
slip-lines and the values of the function ip on the non-characteristic curve OC 
are given. 

Such formulation arises in applied plasticity if a tangential stress r n is 
given on a contact line OC. It is assumed here that r n = /ik, // G [0, 1], 
(Prandtl friction law). Note that constancy of r n implies that the angle r\ 
also remains constant on OC. Assuming that the curve OC is defined by an 
equation f3 = /3(a) = c[a], one obtains elementary identities 

R(a, /3(a))/ cost/ = -S(a, /3(a)) 0/ sinr/ = r(a + (3(a) + r/)(l + (3'). (28) 

The unknown vector c and vector b for the initial /3-line of the equivalent 
initial characteristic problem can be obtained from these identities. A par- 
ticular case is analyzed below. 

Let OC be a straight line. If x-axis is directed along the a-line passing 
through the origin, then yj = a-|-/3 = 0on OC, or equivalently, (3 = —a. 
The rightmost term in (1281) takes the form oo ■ and becomes useless. The 
first identity takes the form 

R(a, — a)/ cosr] = — S(a, — a))/ sinr? (29) 

Using ffl8l) . one obtains after some mathematical manipulations the missing 
boundary data for the initial /3-line of the equivalent initial characteristic 
problem. Specifically, if rj < 7r/2, the row of the coefficients for the /3-line 
radius of curvature is given by the equation b = aT(i]) , where the elements 
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Figure 4: Slip-line field of extrusion/drawing through the short wedge-shaped 
die. 



of the upper triangular matrix T are defined by the equations 



'ii 
ij 




(tan 2 ?] — 1) for j > i, 



It can be shown that for perfectly rough boundary (rj = vr/2) b n = (— l) n a n . 

4 Test problems 

In this section SLAT is applied to two test problems. The first example is 
classical and its exact solution is well known [9]. The second problem was 
treated before only numerically. Since the main mathematical operation used 
in SLAT is multiplication of matrices and vectors of small size, computer time 
needed for each problem is negligible. However, some preliminary analytic 
work is needed to present initial data in the required form. 

4.1 Extrusion through a short wedge-shaped die 

Consider an extrusion of a plastic material through a short wedge-shaped die 
of angle 7 (Fig. 4). For comparability of this example with solution given 
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in [5], the die is assumed to be frictionless (i.e., r\ = 7r/4), and the following 
values of parameters are used: 7 = 10°, ot\ = 30°, /3i = —20°, OA = 2 
(dimensionless units). Thus the initial data are a = (y/2, 0, 0, ...) and b = 
(—y/2, 0, 0, ...). The initial characteristic problem with the data was solved 
in the characteristic rectangle (Fig. 4) using (j!4p . Then the coordinates of 
the point E relative to the origin O were found using simple geometric con- 
siderations. The hydrostatic stress components ps = —a(B), pr> = —a(D), 
and the extrusion pressure P/H were obtained using the solution for the 
initial characteristic problem and Hencky's first theorem in a standard static 
manner [5]. The value of H in this example is found to be 2.28774. 

Table 1. Comparison of some parameters for the extrusion process 





x E 


Ye 


P/H 


[5J 


0.9065 


-2.2877 


0.4117 


SLAT 


0.90648 


-2.28774 


0.41164 



For accuracy of 10 -5 , five-dimensional truncation of matrices and vectors 
was found to be sufficient. The results for some parameters obtained by 
the original matrix method [5] and SLAT are practically identical (Table 1). 
Hydrostatic and extrusion pressures are normalized by the shear yield stress 
of the material. 

4.2 Stress state calculation near loaded elliptic hole 

Consider planar plastic flow near elliptic hole loaded with constant normal 
pressure (Fig. 5). Firstly, Cauchy problem in the domain 1 was solved. 
Secondly, initial characteristic problem was solved in the domain 2 (note that 
the domain 2 is the upper half of a characteristic rectangle). Equation for the 
radius of curvature of the elliptic contour is obtained as follows. Let lengths 
be normalized by the semimajor axis of the ellipse. Then using parametric 
equations x = cost, y = frsint, < b < 1, of the ellipse and the standard 
formula for the curvature, one obtains the radius of curvature pit) in the 
form 

pit) = 2y/2{{\ + b 2 ) + (1 - b 2 ) cos 2tf' 2 /b. (30) 
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Figure 5: Slip-line field for one quarter of an elliptic hole loaded with constant 
normal pressure. 

From elementary geometry it follows that 

7 = 2a + 7r/4 = tan~ 1 (6coU)-7r/4. (31) 

From this equation the variable t can be expressed in terms of a, and after 
some mathematical manipulations one obtains the radius of the curvature 
for the elliptic hole as a function of a : 

r(a) = p(t(a)) = 2^6 2 (1 + d 2 )~ 3/2 (l + q cos4a)- 3/2 , (32) 

where q = (1 — b 2 )/(l + b 2 ). For comparability with [11], value b = 0.4 have 
been chosen. For the domain 1, the row of the coefficients for the initial 
a-line of the equivalent initial characteristic problem was found using the 
first equation in ([25]) . For accuracy of 10 -4 , fifteen-dimensional truncation 
of vectors was found to be sufficient. In general, numerical experimentations 
show that the number of terms needed for a given accuracy grows significantly 
when the value of b decreases. Coefficients for the initial /3-line, which is 
symmetric to the a-line with respect to the horizontal axis, are given by the 
equation bk = (— 1) +1 afc. The main stresses on the £-axis normalized by 2k 
are given by the equations o\ = o + (2a + 1), cr 2 = cr + 2a . The comparative 
results are shown in Fig. 6, where A = a\ + p. 
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Figure 6: Comparison of the computation of A = o\ + p for the elliptic hole 
done numerically [TT] and by SLAT. 

5 Conclusions 

Mathematically accurate and efficient analytic implementation of the slip- line 
method aimed at computing the stress fields for plane strain deformation of 
the rigid-plastic medium has been presented. It can be used as a source of 
reliable test problems for numerical methods and in engineering applications. 

The velocity field in the hodograph plane can also be found by this 
method. Using solutions in the physical and hodograph planes, exact ex- 
pressions for the plastic power can also be derived. This gives an alternative 
way for computing the technological pressure, an important parameter for 
engineering applications, and also provides a framework for solving relevant 
optimization problems. Detailed explicit formulas for computation of ener- 
getic characteristics of slip-line fields will be given elsewhere. 

Although only Prandtl friction law was analyzed in the mixed problem 
presented in this work, Coulomb friction requires just technical modifications 
and can be treated in a similar way. This may be of interest in the theory 
of granular flow, where using this kind of friction law is customary [3J. This 
modification of SLAT will be presented in a later paper. 



1.75 
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